Björn Reineking, 2024-06-19
We will use the buffalo data set.
library(ggplot2)
library(animove)
library(survival)
library(MASS)
library(lubridate)
library(dplyr)
library(nlme)
library(pbs)
library(circular)
library(CircStats)
library(amt)
library(ctmm)
library(terra)
library(sf)
data(buffalo_utm_mv2) # Load buffalo data from animove package
move2_TO_track.xyt <- function(mv2){
if(mt_is_move2(mv2)){
warning("!!INFO!!: only coordinates, timestamps and track IDs are retained")
track(
x=sf::st_coordinates(mv2)[,1],
y=sf::st_coordinates(mv2)[,2],
t=mt_time(mv2),
id=mt_track_id(mv2),
crs = sf::st_crs(mv2)
)
}
}
buffalo_tracks <- move2_TO_track.xyt(buffalo_utm_mv2)
## Warning in move2_TO_track.xyt(buffalo_utm_mv2): !!INFO!!: only coordinates,
## timestamps and track IDs are retained
data(buffalo_env)
buffalo_env <- rast(buffalo_env) # conversion from raster to terra objects
summary(buffalo_tracks)
## x_ y_ t_
## Min. :356761 Min. :7216465 Min. :2005-02-17 05:05:00.00
## 1st Qu.:375417 1st Qu.:7234241 1st Qu.:2005-08-22 06:36:15.00
## Median :379045 Median :7290006 Median :2005-10-14 02:38:30.00
## Mean :380564 Mean :7275531 Mean :2005-10-21 03:03:24.74
## 3rd Qu.:386182 3rd Qu.:7310209 3rd Qu.:2005-12-14 18:32:15.00
## Max. :397937 Max. :7338708 Max. :2006-12-31 14:34:00.00
## id
## Cilla :3528
## Gabs :3670
## Mvubu :2572
## Pepper:4571
## Queen :4545
## Toni :5776
We can plot the buffalo tracks, but they do not have a plot method, so we need to give a bit more information to the points() function.
plot(buffalo_env[[1]])
points(y_ ~ x_, data = buffalo_tracks, col = factor(buffalo_tracks$id))
To speed up analyses, we will only work with one individual, Cilla
cilla <- filter(buffalo_tracks, id == "Cilla")
hist(step_lengths(cilla))
which(step_lengths(cilla) > 5000)
## [1] 1
The very first step is unusually long; let us plot the first day in red on top of the full trajectory.
plot(cilla, type = "l")
lines(filter(cilla, t_ < min(t_) + days(1)), col = "red")
Let us exclude the full first day
cilla <- filter(cilla, t_ > min(t_) + days(1))
It is a good idea to perform the analysis at several temporal scales, i.e. different step durations.
The initial sampling rate of Cilla is about 1 hour:
summarize_sampling_rate(cilla)
## # A tibble: 1 × 9
## min q1 median mean q3 max sd n unit
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 0.05 1 1 0.998 1 2 0.0516 3503 hour
SSF assumes that the velocities are not autocorrelated. So we cannot simply do the analysis at the highest temporal resolution of the data. We could try to use the downweighting trick that we use for the RSF, but for the SSF, the time between successive steps will also affect the point estimates of the parameters, so the problem is a bit more tricky. As a rule-of-thumb, I would downsample the data such that the autocorrelation in velocities has decayed to something between 2% and 5%. Say you had a sampling of 1 hour, and the SSF should be done at 3 hours given this rule of thumb, then you could still use all data, by constructing 3 step data sets, each starting one hour apart, and give each a weight of 1/3 in the likelihood.
cilla_telemetry <- as_telemetry(cilla)
plot(variogram(cilla_telemetry), xlim = c(0, 10 * 3600))
GUESS <- ctmm.guess(cilla_telemetry, interactive=FALSE)
if (file.exists("cilla_fit.rds")) {
FIT <- readRDS("cilla_fit.rds")
} else {
FIT <- ctmm.select(cilla_telemetry, GUESS)
saveRDS(FIT, "cilla_fit.rds")
}
plot(variogram(cilla_telemetry), xlim = c(0, 10 * 3600))
abline(v = FIT$tau["velocity"] * -log(0.01) / 3600, col = "blue")
abline(v = FIT$tau["velocity"] * -log(0.02) / 3600, col = "red")
abline(v = FIT$tau["velocity"] * -log(0.05) / 3600, col = "yellow")
legend("bottomright", lty = 1, col = c("blue", "red", "yellow"), legend = c("1%", "2%", "5%"), title = "Velocity\nautocorrelation",
bty = "n")
Now we resample to 3 hour intervals, with a tolerance of 15 minutes
step_duration <- 3
cilla <- track_resample(cilla, hours(step_duration), tolerance = minutes(15))
Look at the new sampling rate
summarize_sampling_rate(cilla)
## # A tibble: 1 × 9
## min q1 median mean q3 max sd n unit
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 2.95 3 3 3.00 3 4 0.0308 1165 hour
So there is at least two observations that are more than 3 hours 15 minutes apart, so there should be at least two bursts:
table(cilla$burst_)
##
## 1 2
## 322 844
If there are bursts, we may want to filter bursts with very few locations. For example, to calculate a turning angle, we need at least three locations. So we often will want to filter out bursts with at least 3 observations:
cilla <- filter_min_n_burst(cilla, 3)
Convert locations to steps. We will have fewer rows in the step data frame than in the track data frame because the final position is not a complete step.
ssf_cilla <- steps_by_burst(cilla)
We still have steps without a turning angle (the first step in a burst)
which(is.na(ssf_cilla$ta_))
## [1] 1 322
ssf_cilla <- filter(ssf_cilla, !is.na(ta_))
par(mfrow = c(1, 2))
hist(ssf_cilla$sl_, breaks = 20, main = "",
xlab = "Distance (m)")
hist(ssf_cilla$ta_, main="",breaks = seq(-pi, pi, len=11),
xlab="Relative angle (radians)")
fexp <- fitdistr(ssf_cilla$sl_, "exponential")
fgamma <- amt::fit_distr(ssf_cilla$sl_, "gamma")
par(mfrow = c(1, 1))
hist(ssf_cilla$sl_, breaks = 50, prob = TRUE,
xlim = c(0, 8000), ylim = c(0, 2e-3),
xlab = "Step length (m)", main = "")
plot(function(x) dexp(x, rate = fexp$estimate), add = TRUE, from = 0.1, to = 8000, col = "red")
plot(function(x) dgamma(x, shape = fgamma$params$shape,
scale = fgamma$params$scale), add = TRUE, from = 0.1, to = 8000, col = "blue")
legend("topright", col = c("red", "blue"), lty = 1,
legend = c("exponential", "gamma"), bty = "n")
fvmises <- fit_distr(ssf_cilla$ta_, "vonmises")
par(mfrow = c(1, 1))
hist(ssf_cilla$ta_, breaks = 50, prob = TRUE,
xlim = c(-pi, pi),
xlab = "Turning angles (rad)", main = "")
plot(function(x) dvonmises(x, mu = 0, kappa = fvmises$params$kappa), add = TRUE, from = -pi, to = pi, col = "red")
## Warning in as.circular(x): an object is coerced to the class 'circular' using default value for the following components:
## type: 'angles'
## units: 'radians'
## template: 'none'
## modulo: 'asis'
## zero: 0
## rotation: 'counter'
## conversion.circularxradians0counter
## Warning in as.circular(x): an object is coerced to the class 'circular' using default value for the following components:
## type: 'angles'
## units: 'radians'
## template: 'none'
## modulo: 'asis'
## zero: 0
## rotation: 'counter'
## conversion.circularmuradians0counter
Create random steps. We typically get a warning that “Step-lengths or turning angles contained NA, which were removed”, because of the missing turning angles at the start of a burst.
set.seed(2)
ssf_cilla <- steps_by_burst(cilla)
ssf_cilla <- random_steps(ssf_cilla, n_control = 200)
my_step_id <- 4
ggplot(data = filter(ssf_cilla, step_id_ == my_step_id | (step_id_ %in% c(my_step_id - 1, my_step_id - 2) & case_ == 1)),
aes(x = x2_, y = y2_)) + geom_point(aes(color = factor(step_id_))) + geom_point(data = filter(ssf_cilla, step_id_ %in% c(my_step_id, my_step_id - 1, my_step_id - 2) & case_ == 1), aes(x = x2_, y = y2_, color = factor(step_id_), size = 2))
I recommend to always use the option “both”, which provides the environmental conditions at the start and the end of the step. The condition at the end are what we use for selection, and the conditions at the start can be used to modify e.g. turning angles and step length.
ssf_cilla <- amt::extract_covariates(ssf_cilla, buffalo_env, where = "both")
Adding hour modelling diurnal variation in step lengths, turning angles, and preference for environmental conditions
ssf_cilla <- mutate(ssf_cilla, "hour" = hour(t1_) + minute(t1_) / 60)
Remove NA’s
ssf_cilla <- ssf_cilla[complete.cases(ssf_cilla),]
m_1 <- fit_clogit(ssf_cilla, case_ ~ cos(ta_) + sl_ + log(sl_) + slope_end + elev_end +
mean_NDVI_end + var_NDVI_end + strata(step_id_))
summary(m_1)
## Call:
## coxph(formula = Surv(rep(1, 233562L), case_) ~ cos(ta_) + sl_ +
## log(sl_) + slope_end + elev_end + mean_NDVI_end + var_NDVI_end +
## strata(step_id_), data = data, method = "exact")
##
## n= 233562, number of events= 1162
##
## coef exp(coef) se(coef) z Pr(>|z|)
## cos(ta_) 1.029e-02 1.010e+00 4.694e-02 0.219 0.82654
## sl_ -6.329e-06 1.000e+00 6.223e-05 -0.102 0.91900
## log(sl_) 1.632e-02 1.016e+00 3.375e-02 0.484 0.62871
## slope_end -1.113e-02 9.889e-01 2.170e-02 -0.513 0.60786
## elev_end -1.575e-02 9.844e-01 3.241e-03 -4.860 1.18e-06 ***
## mean_NDVI_end 5.743e+00 3.120e+02 2.142e+00 2.681 0.00735 **
## var_NDVI_end -1.391e+00 2.487e-01 9.889e+00 -0.141 0.88811
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## cos(ta_) 1.0103 0.989767 9.215e-01 1.108e+00
## sl_ 1.0000 1.000006 9.999e-01 1.000e+00
## log(sl_) 1.0165 0.983811 9.514e-01 1.086e+00
## slope_end 0.9889 1.011196 9.478e-01 1.032e+00
## elev_end 0.9844 1.015875 9.781e-01 9.906e-01
## mean_NDVI_end 312.0244 0.003205 4.683e+00 2.079e+04
## var_NDVI_end 0.2487 4.020162 9.514e-10 6.503e+07
##
## Concordance= 0.548 (se = 0.008 )
## Likelihood ratio test= 65.46 on 7 df, p=1e-11
## Wald test = 56.6 on 7 df, p=7e-10
## Score (logrank) test = 56.87 on 7 df, p=6e-10
In statistics, multicollinearity (also collinearity) is a phenomenon in which one predictor variables in a multiple regression model can be linearly predicted from the others with a substantial degree of accuracy. In this situation the coefficient estimates of the multiple regression may change erratically in response to small changes in the model or the data.” Wikipedia, accessed 29.08.2017
One way of dealing with collinearity is to select a subset of variables that is sufficiently uncorrelated Dormann et al. 2013. Here we simply look at pairwise correlation between predictors.
round(cor(ssf_cilla[, c("slope_end", "elev_end",
"water_dist_end", "mean_NDVI_end", "var_NDVI_end")]), 2)
## slope_end elev_end water_dist_end mean_NDVI_end var_NDVI_end
## slope_end 1.00 0.37 0.09 0.21 0.10
## elev_end 0.37 1.00 0.77 -0.16 0.58
## water_dist_end 0.09 0.77 1.00 -0.44 0.67
## mean_NDVI_end 0.21 -0.16 -0.44 1.00 -0.40
## var_NDVI_end 0.10 0.58 0.67 -0.40 1.00
elev and water_dist are positively correlated > 0.7
m1_water <- fit_clogit(ssf_cilla, case_ ~ cos(ta_) + sl_ + log(sl_) + water_dist_end + strata(step_id_))
m1_elev <- fit_clogit(ssf_cilla, case_ ~ cos(ta_) + sl_ + log(sl_) + elev_end + strata(step_id_))
AIC(m1_water$model)
## [1] 12332.07
AIC(m1_elev$model)
## [1] 12274.84
So we pick elev, because it by itself explains the movement better.
Fit step selection function
m_1 <- fit_clogit(ssf_cilla, case_ ~ cos(ta_) + sl_ + log(sl_) + slope_end + elev_end +
mean_NDVI_end + var_NDVI_end + strata(step_id_))
summary(m_1)
## Call:
## coxph(formula = Surv(rep(1, 233562L), case_) ~ cos(ta_) + sl_ +
## log(sl_) + slope_end + elev_end + mean_NDVI_end + var_NDVI_end +
## strata(step_id_), data = data, method = "exact")
##
## n= 233562, number of events= 1162
##
## coef exp(coef) se(coef) z Pr(>|z|)
## cos(ta_) 1.029e-02 1.010e+00 4.694e-02 0.219 0.82654
## sl_ -6.329e-06 1.000e+00 6.223e-05 -0.102 0.91900
## log(sl_) 1.632e-02 1.016e+00 3.375e-02 0.484 0.62871
## slope_end -1.113e-02 9.889e-01 2.170e-02 -0.513 0.60786
## elev_end -1.575e-02 9.844e-01 3.241e-03 -4.860 1.18e-06 ***
## mean_NDVI_end 5.743e+00 3.120e+02 2.142e+00 2.681 0.00735 **
## var_NDVI_end -1.391e+00 2.487e-01 9.889e+00 -0.141 0.88811
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## cos(ta_) 1.0103 0.989767 9.215e-01 1.108e+00
## sl_ 1.0000 1.000006 9.999e-01 1.000e+00
## log(sl_) 1.0165 0.983811 9.514e-01 1.086e+00
## slope_end 0.9889 1.011196 9.478e-01 1.032e+00
## elev_end 0.9844 1.015875 9.781e-01 9.906e-01
## mean_NDVI_end 312.0244 0.003205 4.683e+00 2.079e+04
## var_NDVI_end 0.2487 4.020162 9.514e-10 6.503e+07
##
## Concordance= 0.548 (se = 0.008 )
## Likelihood ratio test= 65.46 on 7 df, p=1e-11
## Wald test = 56.6 on 7 df, p=7e-10
## Score (logrank) test = 56.87 on 7 df, p=6e-10
slope and var_NDVI do not contribute significantly to the fit.
Model selection is a vast topic. I recommend using only few models with ecological justification, rather than searching for the “best” model in a huge model space. Here we just use stepwise backward selection based on AIC
m_2 <- fit_clogit(ssf_cilla, case_ ~ cos(ta_) + sl_ + log(sl_) + slope_end + elev_end + mean_NDVI_end + strata(step_id_))
AIC(m_1$model)
## [1] 12273.42
AIC(m_2$model)
## [1] 12271.44
summary(m_2)
## Call:
## coxph(formula = Surv(rep(1, 233562L), case_) ~ cos(ta_) + sl_ +
## log(sl_) + slope_end + elev_end + mean_NDVI_end + strata(step_id_),
## data = data, method = "exact")
##
## n= 233562, number of events= 1162
##
## coef exp(coef) se(coef) z Pr(>|z|)
## cos(ta_) 1.049e-02 1.011e+00 4.691e-02 0.224 0.82307
## sl_ -6.477e-06 1.000e+00 6.222e-05 -0.104 0.91710
## log(sl_) 1.635e-02 1.016e+00 3.375e-02 0.484 0.62812
## slope_end -1.068e-02 9.894e-01 2.147e-02 -0.498 0.61881
## elev_end -1.593e-02 9.842e-01 2.982e-03 -5.342 9.19e-08 ***
## mean_NDVI_end 5.761e+00 3.175e+02 2.140e+00 2.692 0.00711 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## cos(ta_) 1.0105 0.989565 0.9218 1.108
## sl_ 1.0000 1.000006 0.9999 1.000
## log(sl_) 1.0165 0.983784 0.9514 1.086
## slope_end 0.9894 1.010738 0.9486 1.032
## elev_end 0.9842 1.016059 0.9785 0.990
## mean_NDVI_end 317.5217 0.003149 4.7872 21060.322
##
## Concordance= 0.548 (se = 0.008 )
## Likelihood ratio test= 65.44 on 6 df, p=4e-12
## Wald test = 56.41 on 6 df, p=2e-10
## Score (logrank) test = 56.18 on 6 df, p=3e-10
m_3 <- fit_clogit(ssf_cilla, case_ ~ cos(ta_) + sl_ + log(sl_) + elev_end + mean_NDVI_end + strata(step_id_))
AIC(m_3$model)
## [1] 12269.69
summary(m_3)
## Call:
## coxph(formula = Surv(rep(1, 233562L), case_) ~ cos(ta_) + sl_ +
## log(sl_) + elev_end + mean_NDVI_end + strata(step_id_), data = data,
## method = "exact")
##
## n= 233562, number of events= 1162
##
## coef exp(coef) se(coef) z Pr(>|z|)
## cos(ta_) 1.025e-02 1.010e+00 4.691e-02 0.218 0.82707
## sl_ -7.001e-06 1.000e+00 6.223e-05 -0.112 0.91043
## log(sl_) 1.640e-02 1.017e+00 3.376e-02 0.486 0.62706
## elev_end -1.668e-02 9.835e-01 2.592e-03 -6.436 1.23e-10 ***
## mean_NDVI_end 5.390e+00 2.192e+02 2.009e+00 2.683 0.00729 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## cos(ta_) 1.0103 0.989805 0.9216 1.108e+00
## sl_ 1.0000 1.000007 0.9999 1.000e+00
## log(sl_) 1.0165 0.983733 0.9515 1.086e+00
## elev_end 0.9835 1.016820 0.9785 9.885e-01
## mean_NDVI_end 219.2224 0.004562 4.2767 1.124e+04
##
## Concordance= 0.547 (se = 0.008 )
## Likelihood ratio test= 65.19 on 5 df, p=1e-12
## Wald test = 56.28 on 5 df, p=7e-11
## Score (logrank) test = 55.57 on 5 df, p=1e-10
Forester et al. 2009 Ecology 90:3554–3565. Calculate deviance residuals for each stratum (i.e., the sum of the residuals for the case and all associated controls).
ssf_residuals <- function(m, data) {
df <- tibble::as_tibble(data.frame("time" = data$t1_,
"residuals" = residuals(m$model, type = "deviance")))
df <- df |> dplyr::group_by(time) |> dplyr::summarise(residuals = sum(residuals))
df$group <- 1
df
}
resid_df <- ssf_residuals(m_3, ssf_cilla)
Fit an intercept-only mixed-effects model using lme() from the nlme package.
rm1 <- lme(residuals ~ 1, random = ~ 1 | group,
data = resid_df)
plot(ACF(rm1, maxLag = 40), alpha = 0.05)
So we see that there is some significant autocorrelation at lags up to 8.
One effect of residual temporal autocorrelation is too extreme p-values, but it may also cause bias in parameter estimates. Forester et al. 2009 suggest a way to estimate more appropriate confidence intervals and p-values.
Caveat: This habitat preference in general will not match the range distribution of the animal (i.e. how much time it spends where). For the UD, a generic approach is to do simulations from the fitted model (see below). See for more information Signer et al. (2024) Simulating animal space use from fitted integrated Step-Selection Functions (iSSF)
The habitat_preference function assumes that all environmental layers are represented in one raster stack. The habitat_preference function does not work for the “barrier” model.
habitat_preference <- function(model, object, include_avail = TRUE, data_crs = crs(object)) {
remove_end <- function(x) {
names(x) <- gsub("_end", "", names(x))
x
}
beta <- coef(model)
exclude_terms <- c("strata","ta_", "sl_")
for(i in exclude_terms) {
if (length(grep(i, names(beta)))) {
beta <- beta[-grep(i, names(beta))]
}
}
beta <- remove_end(beta)
if(include_avail) {
xy <- as.data.frame(xyFromCell(object, 1:ncell(object)))
xy <- sf::st_as_sf(xy, coords = c("x", "y"), crs = crs(object))
xy <- st_transform(xy, data_crs)
xy <- st_coordinates(xy)
colnames(xy) <- c("x2_", "y2_")
} else {
xy <- cbind("x2_" = rep(0, ncell(object)), "y2_" = rep(0, ncell(object)))
}
newdata <- as.data.frame(cbind(xy, values(object)))
f <- as.formula(paste("~", paste(names(beta), collapse = " + ")))
attr(newdata, "na.action") <- "na.pass"
X <- stats::model.matrix.default(f, data = newdata, na.action = na.pass)
lambda <- X[,-1] %*% beta
r <- object[[1]]
r[] <- exp(lambda - max(lambda, na.rm = TRUE))
r <- r / sum(values(r), na.rm = TRUE)
r
}
habitat_map <- habitat_preference(m_1, buffalo_env)
plot(habitat_map)
lines(ssf_cilla[, c("x1_", "y1_")])
Now zoom in
plot(crop(habitat_map, ext(as_sf(cilla)) + 5000))
lines(ssf_cilla[, c("x1_", "y1_")])
The model is strongly driven by elevation
plot(crop(buffalo_env[["elev"]], ext(as_sf(cilla)) + 5000))
lines(ssf_cilla[, c("x1_", "y1_")])
Here: a model with time-varying preference for mean_NDVI We group observation in 3 hour bins to smooth the picture
boxplot(mean_NDVI_end ~ I(floor(hour/3)*3), data =
filter(ssf_cilla, case_ == 1),xlab = "Time of day",
ylab = "mean NDVI")
We can do the same for other variables, including those of the “movement kernel”, e.g. distance
boxplot(sl_ ~ floor(hour), data =
filter(ssf_cilla, case_ == 1), xlab = "Time of day",
ylab = "dist")
What behavioural rhythm do these figures suggest?
m_time_ndvi <- fit_clogit(ssf_cilla, case_ ~ cos(ta_) + sl_ + log(sl_) + elev_end + mean_NDVI_end +
mean_NDVI_end:pbs(hour, df = 5, Boundary.knots = c(0,24)) + strata(step_id_))
m_time_dist <- fit_clogit(ssf_cilla, case_ ~ cos(ta_) + sl_ + log(sl_) + elev_end + mean_NDVI_end +
sl_:pbs(hour, df = 5, Boundary.knots = c(0,24)) + strata(step_id_))
extract_stratum <- function(object) {
attr(object$terms, 'special')$strata[1]
}
stratum <- extract_stratum(m_time_ndvi$model)
pred_data_ndvi <- data.frame("step_id_" = stratum, ta_ = 0, sl_ = 1, elev_end = 0, mean_NDVI_end = 1, hour = seq(0, 24, len = 101))
m_time_ndvi_clogit <- clogit(case_ ~ cos(ta_) + sl_ + log(sl_) + elev_end + mean_NDVI_end +
mean_NDVI_end:pbs(hour, df = 5, Boundary.knots = c(0,24)) + strata(step_id_), data = ssf_cilla)
pred_time <- survival:::predict.coxph(m_time_ndvi_clogit, newdata = pred_data_ndvi, se.fit = TRUE)
upper <- pred_time$fit + 1.96 * pred_time$se.fit
lower <- pred_time$fit - 1.96 * pred_time$se.fit
par(mfrow = c(1, 1))
plot(pred_data_ndvi$hour, pred_time$fit, type = "l",
ylim = range(c(upper, lower)), xlab = "Time of day",
ylab = "Preference mean_NDVI")
lines(pred_data_ndvi$hour, upper, lty = 2)
lines(pred_data_ndvi$hour, lower, lty = 2)
abline(h = 0, lty = 3)
When using importance sampling in SSF (the default in amt), the parameter estimates for the movement kernel (e.g. step length, log(step length), cos(turn angle)) represent differences to the parameter values of the proposal distribution. Parameter values for the step length and turn angle distributions have to be calculated from the fitted ssf model and the proposal distribution. See Avgar et al. (2016), Appendix 3 for more information.
# Proposal distribution - step length
# sl_distr(m_time_ndvi_fit)
initial_sl_distr <- sl_distr(ssf_cilla)
updated_sl_distr <- update_sl_distr(m_time_ndvi, beta_log_sl = "log(sl_)")
plot_sl_distr <- function(x, n = 1000, upper_quantile = 0.99, plot = TRUE) {
# based on amt:::plot_sl_base
stopifnot(x$name == "gamma")
xx <- x$params
to <- qgamma(upper_quantile, shape = xx$shape, scale = xx$scale)
xs <- seq(0, to, length.out = n)
ys <- dgamma(xs, shape = xx$shape, scale = xx$scale)
if (plot) {
plot(xs, ys, type = "l", ylab = "Probability", xlab = "Distance")
}
invisible(data.frame(sl = xs, d = ys))
}
initial_df <- plot_sl_distr(initial_sl_distr)
updated_df <- plot_sl_distr(updated_sl_distr)
plot(initial_df, type = "l", xlab = "Distance", ylab = "Density")
lines(updated_df, col = "red")
Proposal distribution - turning angle
initial_ta_distr <- ta_distr(ssf_cilla)
updated_ta_distr <- update_ta_distr(m_time_ndvi, beta_cos_ta = "cos(ta_)")
initial_ta_distr$params$kappa
## [1] 0.8178848
updated_ta_distr$params$kappa
## [1] 0.8264717
ta_values <- seq(-pi, pi, len = 101)
plot(ta_values, dvonmises(ta_values, 0, initial_ta_distr$params$kappa), type = "l", xlab = "Turning angle (radians)",
ylab = "Density")
## Warning in as.circular(x): an object is coerced to the class 'circular' using default value for the following components:
## type: 'angles'
## units: 'radians'
## template: 'none'
## modulo: 'asis'
## zero: 0
## rotation: 'counter'
## conversion.circularxradians0counter
## Warning in as.circular(x): an object is coerced to the class 'circular' using default value for the following components:
## type: 'angles'
## units: 'radians'
## template: 'none'
## modulo: 'asis'
## zero: 0
## rotation: 'counter'
## conversion.circularmuradians0counter
lines(ta_values, dvonmises(ta_values, 0, updated_ta_distr$params$kappa), col = "red")
## Warning in as.circular(x): an object is coerced to the class 'circular' using default value for the following components:
## type: 'angles'
## units: 'radians'
## template: 'none'
## modulo: 'asis'
## zero: 0
## rotation: 'counter'
## conversion.circularxradians0counter
## Warning in as.circular(x): an object is coerced to the class 'circular' using default value for the following components:
## type: 'angles'
## units: 'radians'
## template: 'none'
## modulo: 'asis'
## zero: 0
## rotation: 'counter'
## conversion.circularmuradians0counter
The calculation of movement parameters becomes more interesting (challenging) when we have interactions between the movement model parameters (e.g. with environmental conditions). See ?update_gamma or ?update_vonmises
set.seed(23)
k1 <- amt::redistribution_kernel(m_3, map = buffalo_env,
start = amt:::make_start.steps_xyt(ssf_cilla[1, ]),
landscape = "continuous", tolerance.outside = 0.2,
as.rast = FALSE,
n.control = 1e3)
s1 <- amt::simulate_path(k1, n.steps = 500)
extent_tracks <- function(x, y) {
df <- data.frame(na.omit(rbind(x[,c("x_", "y_")], y[,c("x_", "y_")])))
terra::ext(c(range(df$x_), range(df$y_)))
}
elev_crop <- crop(buffalo_env[["elev"]], extent_tracks(s1, cilla) + 2000)
plot(elev_crop)
lines(cilla)
lines(s1$x_, s1$y_, col = "red")
m_hr <- fit_clogit(ssf_cilla, case_ ~ cos(ta_) + sl_ + log(sl_) +
elev_end + water_dist_end + x2_ + y2_ + I(x2_^2 + y2_^2) + strata(step_id_))
summary(m_hr)
## Call:
## coxph(formula = Surv(rep(1, 233562L), case_) ~ cos(ta_) + sl_ +
## log(sl_) + elev_end + water_dist_end + x2_ + y2_ + I(x2_^2 +
## y2_^2) + strata(step_id_), data = data, method = "exact")
##
## n= 233562, number of events= 1162
##
## coef exp(coef) se(coef) z Pr(>|z|)
## cos(ta_) 2.317e-02 1.023e+00 4.712e-02 0.492 0.623
## sl_ 2.601e-05 1.000e+00 6.264e-05 0.415 0.678
## log(sl_) 1.206e-02 1.012e+00 3.376e-02 0.357 0.721
## elev_end -2.156e-02 9.787e-01 3.090e-03 -6.977 3.01e-12 ***
## water_dist_end -1.367e-04 9.999e-01 1.338e-04 -1.022 0.307
## x2_ 1.260e-02 1.013e+00 2.566e-03 4.911 9.04e-07 ***
## y2_ 2.327e-01 1.262e+00 4.784e-02 4.864 1.15e-06 ***
## I(x2_^2 + y2_^2) -1.607e-08 1.000e+00 3.307e-09 -4.860 1.17e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## cos(ta_) 1.0234 0.9771 0.9332 1.1224
## sl_ 1.0000 1.0000 0.9999 1.0001
## log(sl_) 1.0121 0.9880 0.9473 1.0814
## elev_end 0.9787 1.0218 0.9728 0.9846
## water_dist_end 0.9999 1.0001 0.9996 1.0001
## x2_ 1.0127 0.9875 1.0076 1.0178
## y2_ 1.2620 0.7924 1.1490 1.3860
## I(x2_^2 + y2_^2) 1.0000 1.0000 1.0000 1.0000
##
## Concordance= 0.553 (se = 0.008 )
## Likelihood ratio test= 90.33 on 8 df, p=4e-16
## Wald test = 72.22 on 8 df, p=2e-12
## Score (logrank) test = 68.7 on 8 df, p=9e-12
k2 <- amt::redistribution_kernel(m_hr, map = buffalo_env, start = amt:::make_start.steps_xyt(ssf_cilla[1, ]),
landscape = "continuous", tolerance.outside = 0.01, as.rast = FALSE,
n.control = 1e3)
set.seed(2)
s2 <- amt::simulate_path(k2, n.steps = 1000)
plot(crop(buffalo_env[["elev"]], extent_tracks(s2, cilla) + 2000))
lines(cilla)
lines(s2$x_, s2$y_, col = "red")
Plot simulated track on the habitat preference map
habitat_map_hr <- habitat_preference(m_hr, buffalo_env)
terra::plot(crop(habitat_map_hr, extent_tracks(s2, cilla) + 2000))
lines(cilla)
lines(s2$x_, s2$y_, col = "red")
water <- buffalo_env[["water_dist"]] > 100
water <- crop(water, ext(water) - 5000)
plot(water)
ww <- patches(water, zeroAsNA = TRUE)
plot(ww)
ww[is.na(ww)] <- 2
names(ww) <- "water_crossed"
buffalo_env_2 <- c(ww, crop(buffalo_env, ww))
set.seed(2)
ssf_cilla <- cilla |> steps_by_burst() |> random_steps(n_control = 1000) |>
extract_covariates(buffalo_env_2, where = "both") |>
mutate(hour = hour(t1_) + minute(t1_) / 60) |>
na.omit()
m_crossing <- fit_clogit(ssf_cilla, case_ ~ cos(ta_) + sl_ + log(sl_) +
elev_end + water_dist_end + x2_ + y2_ + I(x2_^2 + y2_^2) +
I(water_crossed_end != water_crossed_start) + strata(step_id_))
## Warning in coxexact.fit(X, Y, istrat, offset, init, control, weights = weights,
## : Loglik converged before variable 9 ; beta may be infinite.
summary(m_crossing)
## Call:
## coxph(formula = Surv(rep(1, 1160790L), case_) ~ cos(ta_) + sl_ +
## log(sl_) + elev_end + water_dist_end + x2_ + y2_ + I(x2_^2 +
## y2_^2) + I(water_crossed_end != water_crossed_start) + strata(step_id_),
## data = data, method = "exact")
##
## n= 1160790, number of events= 1160
##
## coef exp(coef)
## cos(ta_) 4.086e-02 1.042e+00
## sl_ 5.968e-05 1.000e+00
## log(sl_) 1.381e-02 1.014e+00
## elev_end -2.112e-02 9.791e-01
## water_dist_end 3.471e-04 1.000e+00
## x2_ 1.059e-02 1.011e+00
## y2_ 2.006e-01 1.222e+00
## I(x2_^2 + y2_^2) -1.388e-08 1.000e+00
## I(water_crossed_end != water_crossed_start)TRUE -1.665e+01 5.848e-08
## se(coef) z Pr(>|z|)
## cos(ta_) 4.711e-02 0.867 0.3858
## sl_ 6.227e-05 0.958 0.3379
## log(sl_) 3.363e-02 0.411 0.6814
## elev_end 3.073e-03 -6.875 6.21e-12 ***
## water_dist_end 1.705e-04 2.036 0.0417 *
## x2_ 2.692e-03 3.932 8.43e-05 ***
## y2_ 4.992e-02 4.017 5.88e-05 ***
## I(x2_^2 + y2_^2) 3.450e-09 -4.022 5.76e-05 ***
## I(water_crossed_end != water_crossed_start)TRUE 6.407e+02 -0.026 0.9793
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95
## cos(ta_) 1.042e+00 9.600e-01 0.9498
## sl_ 1.000e+00 9.999e-01 0.9999
## log(sl_) 1.014e+00 9.863e-01 0.9492
## elev_end 9.791e-01 1.021e+00 0.9732
## water_dist_end 1.000e+00 9.997e-01 1.0000
## x2_ 1.011e+00 9.895e-01 1.0053
## y2_ 1.222e+00 8.183e-01 1.1082
## I(x2_^2 + y2_^2) 1.000e+00 1.000e+00 1.0000
## I(water_crossed_end != water_crossed_start)TRUE 5.848e-08 1.710e+07 0.0000
## upper .95
## cos(ta_) 1.142
## sl_ 1.000
## log(sl_) 1.083
## elev_end 0.985
## water_dist_end 1.001
## x2_ 1.016
## y2_ 1.348
## I(x2_^2 + y2_^2) 1.000
## I(water_crossed_end != water_crossed_start)TRUE Inf
##
## Concordance= 0.564 (se = 0.008 )
## Likelihood ratio test= 140.2 on 9 df, p=<2e-16
## Wald test = 62.23 on 9 df, p=5e-10
## Score (logrank) test = 96.7 on 9 df, p=<2e-16
set.seed(2)
k3 <- amt::redistribution_kernel(m_crossing, map = buffalo_env_2,
start = amt:::make_start.steps_xyt(ssf_cilla[1, ]),
landscape = "continuous", tolerance.outside = 0.01,
as.rast = FALSE, n.control = 1e3)
s3 <- amt::simulate_path(k3, n.steps = 1000)
plot(crop(buffalo_env[["elev"]], extent_tracks(s3, cilla) + 2000))
lines(cilla)
lines(s3$x_, s3$y_, col = "red")
There is a fast method if we have a symmetric and temporally stable jump kernel, e.g. exponential, and no effect of step angles: Barnett, A. & Moorcroft, P. (2008) Analytic steady-state space use patterns and rapid computations in mechanistic home range analysis. Journal of Mathematical Biology, 57, 139–159. The generic but computationally expensive method is to do simulations: Signer et al. (2023) Simulating animal space use from fitted integrated Step-Selection Functions (iSSF). Methods Ecol Evol 15: 43-50
do_steady_state_UD <- TRUE
if(do_steady_state_UD) {
set.seed(23)
s_looong <- amt::simulate_path(k2, n.steps = 1e4) # We should use many more steps here.
xx <- amt::make_track(s_looong, x_, y_, t_, crs = crs(buffalo_env))
xx <- amt::as_telemetry(xx)
ctmm::projection(xx) <- crs(buffalo_env, proj = TRUE)
xx_guess <- ctmm.guess(xx, ctmm(isotropic = TRUE, range = TRUE), interactive = FALSE)
xx_fit <- ctmm.fit(xx, CTMM = xx_guess)
xx_akde <- akde(xx, xx_fit, grid = buffalo_env[[1]])
xx_r <- ctmm::raster(xx_akde, DF = "PMF")
plot(xx_r, col = rev(gray.colors(51, start = 0, end = 1, alpha = NULL)))
}
## Default grid size of 2.98333333333333 hours chosen for bandwidth(...,fast=TRUE).
step_durations <- 3:12
do_run <- TRUE
buffalo_ids <- levels(factor(buffalo_tracks$id))
if(do_run) {
step_interval_simulation_amt <- lapply(buffalo_ids, function(animal) {
lapply(step_durations, function(step_duration) {
ssf_animal <- filter(buffalo_tracks, id == animal) |>
track_resample(hours(step_duration), tolerance = minutes(15)) |>
filter_min_n_burst(3) |>
steps_by_burst() |>
random_steps(n = 200) |>
extract_covariates(buffalo_env) |>
na.omit() |>
filter(sl_ > 0)
m_1 <- clogit(case_ ~ cos(ta_) + sl_ + log(sl_) + elev +
mean_NDVI + strata(step_id_), data = ssf_animal)
list("coef" = coef(m_1), "confint" = confint(m_1))
})
})
}
## Warning in
## random_steps.bursted_steps_xyt(steps_by_burst(filter_min_n_burst(track_resample(filter(buffalo_tracks,
## : Some bursts contain < 3 steps and will be removed
## Warning in
## random_steps.bursted_steps_xyt(steps_by_burst(filter_min_n_burst(track_resample(filter(buffalo_tracks,
## : Some bursts contain < 3 steps and will be removed
## Warning in
## random_steps.bursted_steps_xyt(steps_by_burst(filter_min_n_burst(track_resample(filter(buffalo_tracks,
## : Some bursts contain < 3 steps and will be removed
## Warning in
## random_steps.bursted_steps_xyt(steps_by_burst(filter_min_n_burst(track_resample(filter(buffalo_tracks,
## : Some bursts contain < 3 steps and will be removed
## Warning in
## random_steps.bursted_steps_xyt(steps_by_burst(filter_min_n_burst(track_resample(filter(buffalo_tracks,
## : Some bursts contain < 3 steps and will be removed
## Warning in
## random_steps.bursted_steps_xyt(steps_by_burst(filter_min_n_burst(track_resample(filter(buffalo_tracks,
## : Some bursts contain < 3 steps and will be removed
## Steps with length 0 are present. This will lead to an error when fitting a gamma distribution. 0 step lengths are replaced with the smallest non zero step length, which is: 2.25240718247496
## Warning in
## random_steps.bursted_steps_xyt(steps_by_burst(filter_min_n_burst(track_resample(filter(buffalo_tracks,
## : Some bursts contain < 3 steps and will be removed
## Warning in
## random_steps.bursted_steps_xyt(steps_by_burst(filter_min_n_burst(track_resample(filter(buffalo_tracks,
## : Some bursts contain < 3 steps and will be removed
## Warning in
## random_steps.bursted_steps_xyt(steps_by_burst(filter_min_n_burst(track_resample(filter(buffalo_tracks,
## : Some bursts contain < 3 steps and will be removed
## Warning in
## random_steps.bursted_steps_xyt(steps_by_burst(filter_min_n_burst(track_resample(filter(buffalo_tracks,
## : Some bursts contain < 3 steps and will be removed
## Warning in
## random_steps.bursted_steps_xyt(steps_by_burst(filter_min_n_burst(track_resample(filter(buffalo_tracks,
## : Some bursts contain < 3 steps and will be removed
## Warning in
## random_steps.bursted_steps_xyt(steps_by_burst(filter_min_n_burst(track_resample(filter(buffalo_tracks,
## : Some bursts contain < 3 steps and will be removed
## Warning in
## random_steps.bursted_steps_xyt(steps_by_burst(filter_min_n_burst(track_resample(filter(buffalo_tracks,
## : Some bursts contain < 3 steps and will be removed
## Warning in
## random_steps.bursted_steps_xyt(steps_by_burst(filter_min_n_burst(track_resample(filter(buffalo_tracks,
## : Some bursts contain < 3 steps and will be removed
## Warning in
## random_steps.bursted_steps_xyt(steps_by_burst(filter_min_n_burst(track_resample(filter(buffalo_tracks,
## : Some bursts contain < 3 steps and will be removed
## Warning in
## random_steps.bursted_steps_xyt(steps_by_burst(filter_min_n_burst(track_resample(filter(buffalo_tracks,
## : Some bursts contain < 3 steps and will be removed
## Steps with length 0 are present. This will lead to an error when fitting a gamma distribution. 0 step lengths are replaced with the smallest non zero step length, which is: 1.12625675366255
## Warning in
## random_steps.bursted_steps_xyt(steps_by_burst(filter_min_n_burst(track_resample(filter(buffalo_tracks,
## : Some bursts contain < 3 steps and will be removed
## Warning in
## random_steps.bursted_steps_xyt(steps_by_burst(filter_min_n_burst(track_resample(filter(buffalo_tracks,
## : Some bursts contain < 3 steps and will be removed
## Warning in
## random_steps.bursted_steps_xyt(steps_by_burst(filter_min_n_burst(track_resample(filter(buffalo_tracks,
## : Some bursts contain < 3 steps and will be removed
## Warning in
## random_steps.bursted_steps_xyt(steps_by_burst(filter_min_n_burst(track_resample(filter(buffalo_tracks,
## : Some bursts contain < 3 steps and will be removed
## Warning in
## random_steps.bursted_steps_xyt(steps_by_burst(filter_min_n_burst(track_resample(filter(buffalo_tracks,
## : Some bursts contain < 3 steps and will be removed
## Warning in
## random_steps.bursted_steps_xyt(steps_by_burst(filter_min_n_burst(track_resample(filter(buffalo_tracks,
## : Some bursts contain < 3 steps and will be removed
## Warning in
## random_steps.bursted_steps_xyt(steps_by_burst(filter_min_n_burst(track_resample(filter(buffalo_tracks,
## : Some bursts contain < 3 steps and will be removed
## Warning in
## random_steps.bursted_steps_xyt(steps_by_burst(filter_min_n_burst(track_resample(filter(buffalo_tracks,
## : Some bursts contain < 3 steps and will be removed
## Warning in
## random_steps.bursted_steps_xyt(steps_by_burst(filter_min_n_burst(track_resample(filter(buffalo_tracks,
## : Some bursts contain < 3 steps and will be removed
## Warning in
## random_steps.bursted_steps_xyt(steps_by_burst(filter_min_n_burst(track_resample(filter(buffalo_tracks,
## : Some bursts contain < 3 steps and will be removed
do_run <- TRUE
if (do_run) {
for (j in seq(step_interval_simulation_amt)) {
model_list <- step_interval_simulation_amt[[j]]
name <- names(split(buffalo_utm))[j]
coefs <- sapply(model_list, function(x) x$coef)
ci_lower <- sapply(model_list, function(x) x$confint[, 1])
ci_upper <- sapply(model_list, function(x) x$confint[, 2])
par(mfrow = c(3, 2), mar=c(4,5,1,1), oma = c(1,1,5,1))
for (i in rownames(coefs)) {
plot(c(0,0), xlim = range(step_durations),
ylim = range(c(0, ci_lower[i, ],
coefs[i,],
ci_upper[i, ])), type = "n", xlab = "Step duration (hr)",
ylab = i)
abline(h = 0, lty = 2, col = "red")
lines(step_durations, ci_lower[i, ], lty = 3)
lines(step_durations, ci_upper[i, ], lty = 3)
lines(step_durations, coefs[i, ])
}
mtext(name, outer = TRUE)
}
}
## Loading required package: move
## Loading required package: geosphere
##
## Attaching package: 'geosphere'
## The following object is masked from 'package:amt':
##
## centroid
## Loading required package: sp
##
## Attaching package: 'sp'
## The following object is masked from 'package:amt':
##
## bbox
## Loading required package: raster
##
## Attaching package: 'raster'
## The following objects are masked from 'package:ctmm':
##
## projection, projection<-
## The following object is masked from 'package:amt':
##
## select
## The following object is masked from 'package:nlme':
##
## getData
## The following object is masked from 'package:dplyr':
##
## select
## The following object is masked from 'package:MASS':
##
## select
##
## Attaching package: 'move'
## The following object is masked from 'package:ctmm':
##
## speed
## The following object is masked from 'package:amt':
##
## speed